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ABSTRACT 


A transient overpower (TOP) accident in a Liquid Metal 
Fast Breeder Reactor (LMFBR) is considered. The analysis 
is formulated to model the dynamic response of the reactor 
fuel subassembly during the initial period of the postulated 
Overpower transient. An equivalent cylindrical cell is used 
to model the fuel subassembly. The governing neutronic and 
heat transport equations for each region (fuel, clad, and 
coolant) of the equivalent cylindrical cell are developed. 
Nuclear Doppler broadening feedback is included in the dynam- 
ic model making the coupled equations non-linear. The re- 
Sulting non-linear partial differential field equations are 
transformed into a system of ordinary differential equations 
by the finite element method. An isoparametric, quadratic, 
rectangular element is used for the discretization of the . 
Spatial domain. When using the finite element method, large 
System matrices may result. To facilitate solution of these 
large systems, an optimum compacting scheme is utilized. 
The implicit Gear's method is used for the solution of the 
System of ordinary differential equations. The results for 


a sample problem are presented. 
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I. .2NTRODUCTION 


As the world's fossil fuel resources are depleted, more 
emphasis is being placed on the breeder reactor as a poten- 
tial means of solving the coming energy crisis. While the 
development of new energy sources is being pushed, equal ef- 
fort is being given to the maintenance of an environmental- 
ly clean world. To this end, the safety of breeder reactors 
is receiving a considerable amount of attention before assum- 
ing that the breeder reactor is the answer to the energy 
problem. 

The Liquid Metal Fast Breeder Reactor (LMFBR) appears to 
be one of the most promising breeder reactors. Most engineers 
will concede there is little probability of a nuclear explo- 
Son occurring in the operation of a nuclear reactor. Of 
major concern to engineers is the loss of coolant accident 
(LOCA) and the transient overpower accident (TOP). The pres- 
ent analysis iS concerned with a TOP accident in a LMFBR. 

The analysis is formulated to model the dynamic response of 
the reactor fuel subassembly during the initial period of the 
postulated overpower transient. The primary consideration is 
given to the early response of this fuel subassembly to 
various conditions of disturbances. The phenomenon which oc- 
eurs after core disassembly (i.e., clad melting) is not the 
concern of this analysis. Only the time prior to clad melting 


is being considered. 
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No consideration is given here as to how the overpower 
transient occurs or to why the safety features of the reac- 
ton did not operate properly. It ts postulated that the 
accident has occurred. In this analysis, the TOP accident 
is created by either a step increase in reactivity, a ramp 
increase in reactivity, or a combination of both. 

An inherent safety feature of most reactors, nuclear 
sepia broadening feedback, is included in the dynamic model 
of the fuel subassembly. The Doppler feedback acts to reduce 
the effect of the excursion. Consideration of this feedback 
creates a non-linear system model which is described by a 
non-linear, initial-boundary-value problem. 

The conventional method of solution uses the standard 
point kinetics formulation. Recent studies have pointed out 
a non-negligible error in this model [1], particularly with 
asymmetric disturbances [2], or space-dependent feedback [3]. 
In Ref. [4], a somewhat novel approach of using the finite 
element method (FEM) for the space-time dependent solution 
of the reactor dynamics problem was demonstrated. The FEM 
is effective in handling these asymmetric disturbances and 
Space-dependent feedbacks. Therefore, the finite element 
method was used so that the spatial effects on the postu- 
lated problem may be studied further. 
| The purpose of this work was to demonstrate further the 
applicability of the FEM to the non-linear reactor dynamics 
Problem as well as to investigate the dynamic response of the 


reactor fuel subassembly. The analysis required a novel 
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approach to handle the gap conductances present at the in- 
terfaces of the equivalent cell model of the fuel subassem- 


bly; this will be clarified in the analysis. 


LS 


Ti, DESCRIPMION OF PROBLEM 


Me PHYSICAL SYSTEM 

The typical Liquid Metal Fast Breeder Reactor (LMFBR) 
core consists of many hexagonal modules, each containing 
several hundred fuel pins. For this analysis, an equivalent 
cylindrical cell is used to model the fuel subassembly; see 
Figure 1. The use of equivalent cells as models for larger 
Systems has been common practice in nuclear analysis (i.e., 
the well known Wigner-Sietz method). In using an equivalent 
cell, the actual shape of the reactor core is not important, 
and the analysis is applicable to any reactor which has the 
Same equivalent cell. 

The equivalent cell considered in this analysis, Figure 1, 
is fueled with enriched uranium dioxide, has a stainless 
steel cladding, and has liquid sodium for a coolant. The 


dimensions used are 


a = 0.254 cm, 

bi 0.292 -cm, 

Cii=| 0.4365 ‘ch, 
and i i= 30.0: Cit. 


The gap between the fuel and cladding is very small and, in 
Fact, may be nonexistent as in bonded fuels. The dimension 
of this gap has been assumed negligible. The height, H, of 
the fuel rod is shorter than many proposed systems (Fast Flux 


Testing Facility and Clinch River Breeder Reactor). However, 
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Figure 1. Equivalent Cylindrical Cell 
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to facilitate the numerical solution a smaller rod was used. 
The dynamic behavior prior to fuel pin failure for this sys- 


tem should be similar to the behavior of larger systems. 


The treatment of this problem in three dimensions would 
be prohibitive in computer usage. Therefore, azimuthal sym- 
metry is assumed, and the problem becomes a two-dimensional 


cylindrical (r,z) problem. 


Be SYSTEM MODEL 

The analysis considers the monoenergetic neutron diffu- 
Sion approximation to model the transient neutron transport 
problem. A simple conduction-convection heat transfer model 
is used for the energy transport problem. The tenperarnene 
in the model are directly coupled to the neutron flux through 
the nuclear heat generation within the fuel. The neutron 
population is in turn coupled to the temperature through any 
of a number of reactivity feedback mechanisms. The nuclear 
Doppler effect is perhaps the most important of these mecha- 
nisms since it provides a negative temperature coefficient 
which increases the inherent stability of the reactor. , Prior 
to core disassembly and fuel melting, the nuclear Doppler 
effect is the most dominant feedback and is, therefore, the 
only feedback mechanism considered in this analysis. 

For irradiated, mixed-oxide fuels, a phenomenon of fuel 
restructuring has been commonly observed. This restructuring, 
essentially a change of phase of the fuel material, presents 
a unique heat transfer problem particularly during transient 


conditions. The problem has not been fully characterized and 
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is beyond the scope of this analysis. The fuel Weer there- 
fore, assumed to be a homogeneous mixture of enriched uranium 
dioxide. 

At the fuel-cladding interface, there exists a gap which 
produces a thermal resistance. This thermal resistance is 
One of the most significant deterrents to the energy transfer 
to the coolant. The interface may be in physical contact or 
an actual gap may exist. The prediction of the thermal resis- 
tance is extremely complicated and must take into considera- 
tion many parameters: initial dimensions, type of bond, fill 
gas composition, fuel restructuring, fuel swelling, prior fuel 
life cycle, to mention a few. Reference [5] documents a com- 
puter program which attempts to predict the gap conductance, 
Noap- This treats the thermal resistance at the interface in 
the same manner as a convection heat transfer coefficient when 
considering convection heat transfer. Since it is not the 
objective of this analysis to predict the gap coefficient, a 
representative set of values for gap coefficient, as given in 
Ref [6], is used in this analysis. These values are assumed 
to remain static during the transient. The gap conductance 
Peiiie actually varies with time and will have an effect upon 
Ene transient, as noted in Ref. [7]. The prediction of this 
variance was not considered important for this analysis; 
therefore, the static assumption was made. 

An average convection heat transfer coefficient was used 
to determine the heat transfer from the cladding to the cool- 
ant. The value used was determined from an empirical formula 


given in Ref. [5] and repeated in Appendix C. 
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In this work, consideration is given to step and ramp in- 
creases in reactivity, although any reactivity transient may 
easily be considered. The step and ramp increases in reactiv- 
ity probably represent the most realistic physical reactivity 
mipuecs in a reactor. Once the reactivity has been inserted, 
the transient overpower excursion begins. Unless the Doppler 
feedback can override the inserted reactivity, the excursion 


will continue until there is physical core disassembly. 


C. NUMERICAL SOLUTION 

The system of equations which models the proposed problem 
is a non-linear, initial-boundary-value problem. The conven- 
tional method of solution of the reactor dynamics is the point 
minepies formulation. It was pointed out in Refs. [1], [2], 
[3], and [4], that there is a non-negligible error in this 
model, particularly under conditions of asymmetric disturbances 
Or space-dependent feedback. Reference [4] demonstrates the 
somewhat novel approach of using the finite element method 
(FEM) to solve the space-time dependent reactor dynamics prob- 
lem. As shown in Ref. [4], the FEM is quite effective in 
handling localized perturbations and space-dependent feedback. 
In this work, only uniform disturbances were considered; how- 
ever, the feedback model was space-dependent. Therefore, the 
finite element method is used to solve the non-linear, coupled, 
Space-time dependent neutronic and heat transport field equa- 
tions. The solution technique results in a large computer 
Storage requirement; therefore, an optimum compact storage 


Scheme, Ref. [8], is utilized for storage of the discretized 


matrices. 
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Once the domain has been discretized by the FEM, the solu- 
tion of the resulting ordinary differential equations was to 


be accomplished by the implicit Gear's method, Ref. [9]. 


PAS 


Pris” MODER  PEVELOPMENT 


A. NEUTRONIC ANALYSIS 

In this section the governing field equations for the 
neutron population (flux) for each of the three regions 
(fuel, clad, and coolant) of the domain will be formulated. 
The monoenergetic, diffusion theory will be used. 

Consider an arbitrary volume of material within a reac- 
tor. Applying the condition of conservation to the mono- 
energetic neutrons leads to the neutron equation of 
Continuity [10] 


=-1S. hs E)) = Z(r)o(r,t) = div JdiGr,t) (1) 


where 


G = Spatial. point 
t - time 
- neutron density 


nies 3) 
S(r,t) - neutron production 


z(r) - neutron absorption cross section 
o(r,t) - neutron flux 
J(r,t) - neutron current 


The left-hand side of equation (1) represents the time rate 
of change of the neutron density which is related to flux, 


eo ;, by 


where 
v - neutron velocity 


On the right-hand side of equation (1), the first term is 
neutron production, the second term is a neutron loss through 
absorption, and the third term is a neutron loss through leak- 
age from the control volume. 

Using equation (2) and applying Fick's Law to the equation 


of continuity results in the classical neutron diffusion equa- 


tion 
Rees D(r) - neutron diffusion coefficient 


The vector notation used here is intended to include only two 
dimensions (r,z) since azimuthal symmetry has been assumed. 

Equation (3) is applicable to each of the three regions 
Giethe equivalent cylindrical cell. The subscripts F (fuel), 
c (cladding), and co (coolant), will be used to denote these 
regions. 

1. Fuel Region 

In applying equation (3) to the fuel region the 

material properties of the fuel must be used. Within the 
fuel the neutron souce term is due to the nuclear fission 
Process. During fission, neutrons are released as both 


prompt neutrons and delayed neutrons so that 


2) eee slip a lhie a alr (4) 
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where S, (rst) - prompt neutron source 


p(t» t) - delayed neutron source 


The neutron sources are commonly represented as [10] 


S, =k, Zap op (1-8) (5) 
and 
n 
Sp Fd Du C54, (6) 
i=] 
where 
kK, ~ infinite multiplication factor 
8B - fraction of fission neutrons which appear 
as delayed neutrons 
wn - number of delayed neutron groups 
i 
C. - concentration of delayed neutron precursors 


in the ith group 


AX. - decay constant of the delayed neutrons 


The space and time variables, r and t, will be dropped except 
where needed for clarification. 

The concentration of delayed neutron precursors, C.s 
is given by the following first order partial differential 


equation [10] 


oC. 
hte 
yy Ss alge Ps (7) 
where 
8B. - fraction of delayed neutrons which appear 


as delayed neutrons in the ith group 
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aril 


ther solution of equation (7) is 


t 

= -A.(t-t') ! ! on=Azt 

C. = Bf Bi k,(rot')Zspopdt'+ Coe "1 (8) 
(@) 


C° - concentration of delayed neutron precursors 
Of the 1th group at time zero 


Reference [10] develops an expression for the initial concen- 


tration of delayed neutrons 


of = 9 ° 
Pt see ae i) 


where 


ko = initial finite multiplication factor 


° - initial neutron flux 


companang equations (3), (4), (5), (6), (8), and (9), yields 


the governing equation for the fuel region 


Ve (DeVoe) 5 Dp oglk(1-8) -1] 


n 
myo By a tytt-e L eouunsialh) eric 
=| 0 
oo 
, jie ee 
Eo eee irae ee 


2. Cladding Region 


The cladding separates the fuel and coolant and con- 
tains the fuel and the fission by-products. Equation (3) 


Governs the neutron flux im the clad. Since the clad contains 


2A) 


no fissile material, the neutron source term is zero. The 


field equation is, then, 


V-(D.Vo.) = 


3. Coolant Region 


The annular region around the cladding of the equiva- 
lent cell contains the coolant. As in the clad, there is no 
fissile material in the coolant and, therefore, no neutron 
source term. The equation governing the neutron flux in the 


coolant is 


2 
co Be! - Laco@eco Ma.vi at (12) 


Bauations: (0 )y0 Cl Wencand (he2)) - are the: one-veloci ty, 
diffusion approximation used to model the neutron transport 
problem. 


4. Infinite Multiplication Factor 


Phe a ntinite anu litip bicataon factor, k may be ex- 


pressed as the infinite multiplication factor at time zero 
fsctart of transient), ket pilus thelpos tu lated greaictavity tn- 
sertion (such as a step or a ramp), 9, minus the change in 
the Doppler reactivity feedback, Adp- Other feedback mechan- 
isms are normally not as significant as the Doppler broaden- 


ing feedback prior to fuel melting and have been neglected 


in this analysis. Therefore, 


co 


eset > 8D, (13:) 
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For a fast reactor, k> may be approximated as 


2 
ko =v ge (14) 
aF 
where 
Vv - average number of neutrons released 


per fission 
Ler = FT1ISSTON Cross section of the fuel 


DAE - absorption cross section of the fuel 


The nuclear Doppler effect is a very important safety 
feature in a nuclear reactor. Nuclei in an atom are in con- 
tinual motion due to their own thermal energy. As a result 
of this motion, even when monoenergetic neutrons interact 
with the atom, there appears to be a spread in the energy of 
the neutron - the Doppler effect. It can be shown that the 
cross section of a resonance becomes less in magnitude and 
wider as the motion of the nuclei increases [10]. As the 
temperature increases, the motion increases, and the shape 
of a resonance cross section broadens. This broadening in- 
creases the average cross section, thus, providing a negative 
temperature coefficient. This effect is shown in Figure 2. 
It is this nuclear Doppler broadening effect which provides 
One of the few inherent, reliable, negative reactivity feed- 
backs which slows an overpower transient and possibly stops 


a mild overpower transient. 
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CROSS SECTION (barns) 


Ergure 2, 


ENERGY [ev] 


Doppler Broadening of a Resonance Peak 


28 


The Doppler reactivity change with respect to fuel 


temperature changes should be written as [6] 


QO 


p * 
D aT B/i2 


d 


-] 


bios cr 


(15) 


where a, b, and c are parameters determined from experi- 
mental work and m is an integer. However, as noted in 
Ref. [6], a substantial amount of work has shown that 

T aT is very nearly constant over the temperature range 
under consideration. Therefore, the coefficients a and 


c have been set equal to zero, and b is defined as 
b= Ky, = 7 = (16) 


The constant, Kp is commonly called the Doppler constant. 


Solving equation (16) for Pp yields 


=i i) Dew tk (17) 


°D F 


where K is an integration which may be obtained from 


migetal conditions. 


where 


Pn 7 Doppler effect at time zero 


Te - fuel temperature at time zero 


XS) 


re 


Substituting for K in equation (17) will give 


One a ep App = b In(T-/TP) (19) 


The infinite multiplication factor now becomes 


Keegan 205 = 5 In(Te/Te) (20) 
The effect of delayed neutrons is small compared to 
the prompt neutron effect; therefore, it may be assumed that 
the Doppler effect on delayed neutrons is insignificant to 
tnesoveral!l problem... The .k, of equation (8) is, then, 
assumed to be k2 . With this assumption and equation (20), 


the neutron diffusion equation in the fuel, equation (10), 


may be rewritten as 


V+(DeVoe) + Fide lk2(1-8)-140(1-B8)-b(1-8) In(Te/T2) J 
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To facilitate the present analysis, the number of 
delayed neutron groups is taken as one averaged group. So 
Enat, r; and Bs become X} and 8, respectively. This approxi- 
mation should have little effect on the problem under con- 
Sideration. 

9. Boundary Conditions 

The neutron diffusion problem involves solution of 

tne partial differential equations (11), (12), and (21), with 


the following boundary, interface, and initial conditions: 
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5) De i cle (biz 7) D a6 ar (h52s,,t ) 
a 
6) z(c,z,t) = 0 
Pree RE PAGS Sees Sih ee tah) A inne een 
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Boundary condition 1) results from the assumed azimuthal sym- 
merry. Interface ‘conditions 2); 3), 4), and 5) are continu- 
ity conditions of the flux. Boundary condition 6) results 
from the use of an equivalent cell and basically indicates 
there is an equal number of neutrons transferred in and out 
of the cell at the outer boundary. This should be valid un- 
less the cell is located near the outer edge of the reactor. 
Boundary condition 7) is an assumption that the flux is zero 
at the axial boundaries of the cells. Jnitial conditions 6). 
oy. and 10) are .the assumed initial distributions of, the 


neutron flux. 
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B. HEAT TRANSFER ANALYSIS 

In this section, the principle of conservation of energy 
Will be used to formulate the governing field equations for 
the heat transport in each of the three regions. A simple 
heat conduction model with convection heat transfer to the 
coolant is used to model the heat transport problem. A gap 
conductance model is used to describe the heat transport 
across the gap at the fuel-clad interface. 

1. Fuel Region 

Conservation of energy within the fuel region yields 


the unsteady heat conduction equation with a generation term 
Ve(ke(r)vTE (r,t) )+4(r,t)=pp(r)C 


where 


k-(r) - thermal conductivity of the fuel 
c(r,t) - fuel temperature 


Q(r,t) - nuclear energy generation per unit 
volume 


pe(r) - fuel density 


(r) - fuel specific heat 


As in the neutronic analysis, the vector notation is intended 
to include only two dimensions, (r,z). The r and t will be 
dropped except where needed for clarification. 


The nuclear generation term may be expressed as 


Q =e Tgp of (23) 
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where e - nuclear energy released per fission 


As can be seen, it is through the nuclear generation term 
that the temperature is directly coupled to the neutron 
flux. This coupling and the temperature dependent Doppler 
reactivity feedback combine to make the coupled problem non- 
linear. 

Substituting equation (23) into equation (22) yields 
the governing thermal equation for the fuel 


aT. 


Cex (24) 


Ve(kevT pF Ot 


pute) 


ter = 


fF °F ” PF 


e. Cladding Region 


Conservation of energy within the clad will yield 
the heat conduction equation. With nuclear generation, a 
relatively small amount (~5%) of the energy will be released 
in the cladding and the coolant. However, in this analysis 
it is assumed that the total energy release is in the fuel 
region. This should not create any Significant error, Us- 
ing this assumption, the unsteady heat condition equation 


for the cladding becomes 


3. Coolant Region 


Once again, conservation of energy will lead to 


the heat conduction equation plus an additional term which 


33 


takes into consideration the coolant flow. The governing 


equation is 


oT 
‘ . ok, = co 
Dkeg ven! = “co 3z "co pco co! Peo “nco BE 
(26) 
where Mises - coolant flow velocity 


Equations (24), (25), and (26) are the governing 

equations used to model the energy transport problem. 
4. Interface Conditions 

The interface between the fuel and the cladding may 
be an actual gap with a finite distance, or the surfaces 
may be in intermittent contact on a microscopic scale. To 
model the heat transfer across this interface, a gap heat 
transfer coefficient is introduced. The gap coefficient 
must take into consideration many items (e.g., radiation 
heat transfer across the gap, heat transfer by solid-to-solid 
contact, heat conduction across a gas filled gap). The pre- 
diction of this gap coefficient is extremely complicated and 
beyond the scope of this work. In Ref. [6], a set of values 
for Hap 1s given and the axial variation of Hoap in thas 
analysis is approximately the same. A cosine curve has been 
fitted to the sample data to determine the gap coefficient, 
see Figure 3. The heat flux across the fuel-clad interface 


is, then, 
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The heat conducted out of the fuel is governed by Fourier's 
equation and is equal to the heat transferred across the gap 


ae 
aa eye (28) 


Equating equations (27) and (28) gives the fuel-clad inter- 


face condition 


and from continuity 


ate al 
kelas2) = -(anz.0) > ki (a5 2) = (a. Z5t) (29a) 


As is common practice in convection heat transfer 
amialysis, a coolant surface heat transfer coefficient, hsurf, 
is used to account for the thermal resistance at the clad- 
coolant interface. Similar to the fuel-clad interface analy- 


sis, the clad-coolant interface conditions may be determined 


ee ay (Bz, 2) ate) 
T.(b,z,t) + ~——— —(b,z,t) = T Dysizaant 30 
G cae or oe 
and 
oT. at 6 
k (b»zZ) sp (b.z,t) = k 69 (b»2Z) or (buiz,t.) (30a) 


5. Boundary Conditions 


The boundary and initial conditions for the heat 


transport problem are: 
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at es 
4) ar be sizet) = 0 
oT oT 
F H = Cc H = 
5) ioe (Gerace zet) = oan (n= 72 t) = i) 
6) Terese) = oTe(r) 
7) Te tisO) = Gi) 
8) T.(ps0)h = ote atm) 


Boundary condition 1) results from the assumed azimuthal sym- 
metry. The coolant has been assumed to enter the flow 
Channel at a constant temperature, ToLENUM®* This results in 
condition 2). Boundary conditions 3) and 5) result from an 
assumption that no heat is transferred axially from the fuel 
rod. Boundary condition 4) is the result of the use of the 
equavalent cell. Conditions 6), 7), and 8) are the assumed 
iwenal conditions. 

Solution of the nonlinear coupled neutronic and 
energy transport problem involves the solution of the partial 
aeererential equations (11), (12), (21), (24), (25), and 
(26), with the appropriate boundary, interface, and initial 
conditions. 

Several works, Refs. [4], [£11], [12], have demon- 
Strated the feasibility and success of the finite element 


method in solving nuclear reactor dynamics problems. The FEM 


is used to reduce the partial differential equations developed 
in this analysis to a system of ordinary differential equa- 
tions. Integration of these ordinary differential equations 


(ODE) yields the solution. 
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IV. FINITE ELEMENT FORMULATION 


In this section, the basic theory underlying the finite 
element method is formulated. Selection of the finite ele- 
ments and the shape functions for the elements are given. 
Some simple transformations which facilitate the integration 


necessary in the FEM are also presented. 


A. BASIC THEORY 

To obtain a numerical solution, the governing partial dif- 
ferential field equations are transformed into a system of 
ODE in finite dimensional vector space. This may be accom- 
plished in several manners such as the finite-difference 
method, the variational method, or the weighted residual 
method. In this work, the Galerkin method (a weighted re- 
Sidual method) is utilized for the discretization of the 
Spatial domain. The Galerkin procedure will be applied to 
each of the governing field equations. Any of these equa- 


tions may be considered to be in the following form 


oU rst) = Pvlest) = (rst) (31) 


where w represents the unknown function, e.g., in equation 
(21), ww represents dr >» f represents the operator for 
eden individual equation,.and f is a. foreing function...in 


the finite element method the solution is approximated as 


N 
p(r,t) = Vie. ©) ee N.(r) vy; (t) = <N.>{p, (32) 


a9 


where N - the number of degrees of freedom 

N. - the element shape function 

v; - unknown coordinate function 

< >= Matrix notation for a row vector 

{ } - matrix notation for a column vector 


<N.> - aN No ae WN oo ENS 


tite residual, R(r,t)., is a measure of the error in this 
finite element approximation. The residual may be considered 


as 


_ ow 
RF ae =u = 5 (33) 


The best solution for % is one which "minimizes" this 
residual. Various "minimums" are obtained by the weighted 


residual method by setting 


fw le) R av = 0 fen (34) 
V W, (1) - weighting functions 

With the Galerkin method, the weighting functions are the 
shape functions defining the approximation of equation (32) 


Cree.) W? = N;). A noteworthy attribute of the Galerkin 


method is the opportunity of using an integration-by-parts 
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of the terms involving the second order spatial derivatives. 
A lower order finite element may be used than would have 
been possible otherwise. Once the weighting functions have 


been chosen, the problem becomes 


e) : 
if N; (Se - gy rr ft dV = 0 NICER. che (35) 
V 


The integration involved in equation (35) is carried out 
On the element level, taking advantage of the use of a “local” 
coordinate system. Once the integration is accomplished, 
the results are merged into a system using "global" coordi- 


nates. On the element level 
iY = len PR ceuie Gees. AN (36) 


where the superscript e indicates the element level. Sub- 
stituting ze into equation (35) and noting {yj} iis) aot 


aeumetion of the spatial domain yields 


- CrP p32. 1 e @ - ee é e _ 
J<Ngrr iN Pave Cd IP -f [<Ny>> LANs} -<Ny> fo ]dve (yj 3° = 0 
V V 

(37) 
where Tiise ieee NG 
N° - number of degrees of freedom for 
an element 
The operator L will vary depending upon which governing 


equation is under consideration. 
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Bee SHAPE FUNCTIONS 


The shape functions, N are chosen to satisfy certain 


as 
completeness and convergence criteria [13] and will depend 
upon the finite-element used for the spatial discretization. 
Many previous works, Refs. [4], [11]. and [12], utilized 
linear triangular shaped elements to discretize the spatial 
domain. This element was the first element considered. 
However, because the width of the cladding is very thin, 
elements in the cladding region would have extremely large 
aspect ratios (ratio of base to height) unless an extremely 
large number of elements in the axial direction were used. 
A large number of elements becomes numerically untractable. 
Previous experience with triangular elements had shown that 
large aspect ratios yield inaccurate results. Zlamal, Ref. 
[14], showed the error, e, when using triangular elements, 
i proportional to the square of the longest side, h, and 


inversely proportional to the sine of the smallest angle, y 


erwin yen. 


A triangular element with a large aspect ratio necessarily 
must have a small related angle which adversely affects the 
error in the FEM. Hopefully to alleviate the problem, an 

isoparametric, quadratic, rectangular element was selected. 
The aspect ratio would still be large, but experience, Ref. 
[15], with the use of rectangular elements indicated that a 


large aspect ratio is not always a detrimental factor. 
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The shape functions for this element are well documented, 
Ref. [13]. Utilizing a local" coordinate system (See Figure 


4), the shape functions may be written as 


Corner nodes TS? 135597 
Ne = g(14E,)(14n,) (E5tng-1) (38) 
Midside nodes N, = 5(1-€*)(1#n,) i=2,6 (38a) 
N, = g(1tE,)(1-n?) i=4,8 (38b) 
where Spe Ee 
No = nn; 


These normalized shape functions are shown in Figure 5. 
The local and global coordinates are related by the 


following 


and z= <N3>° a (39a) 


C. COORDINATE TRANSFORMATIONS 
When using a local coordinate system, some simple trans- 
formations facilitate the integrations required by equation 


(37). In cylindrical coordinates with azimuthal symmetry 
dV = 2nrdrdz (40) 


The derivative terms may be transformed by the following 


oN. ON. 

pT Heol 

or ee 

ant? EY] Jon, ay 
Bical weed 

dZ an 
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Figure 4 .v 9 Element Transformation 
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Figure 5. Normalized Shape Functions 
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where aan. is the inverse of the 2x2 Jacobian matrix de- 
fined in Appendix A. As shown in Appendix A, this inverse 


can be easily shown to be (for this problem) 


x * 
a 1 
-] 
[J] = (42) 
* * 
Jo, 922 
where 
* an 
eh = aA IN a Yio = 0 
ee a a 
Joy = 0 P Joo =e) Ne 
AY - radial length of the element 
AZ - axial length of the element 
Elements of area transform as 
drdz = det[J] d&dn (43) 


For this particular problem, det[{J] may be shown to be 
(Appendix A) 


det[J] = = (44) 


where Ac - area of the element 


Elements of axial length become 
e 
pay 


where bo =) axial length of the element 
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Utilization of these transformations makes the integrations 
required in equation (37) amenable to integration by numeri- 


cal Gaussian quadrature. 
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V. APPLICATION OF FEM TO GOVERNING 
FIELD EQUATIONS 
In this chapter, equation (37) is applied to the govern- 
ing field equations previously derived. The element matrices 
for each of the operators are developed so that the discreti- 
zation of the spatial domain may be accomplished. The inte- 
grations required by the application of equation (37) are 


performed numerically using Gaussian quadrature. 


A. GAUSSIAN QUADRATURE 
Prior to the application of equation (37), it is appropri- 
ate to discuss briefly the procedure used for the numerical 


integration. The product Gaussian quadrature formula is [16] 


Te m m 
b= Wt pcenlende = oe WW, G(Eion,) (46) 
As Tea 12) 921 
where Ty - area integration 


g(@jn) = any function of GE vand 7 


- weight associated with location 
tor Jj 


m - number of Gauss sampling points 
in one-dimension 
The values of the weights associated with each Gauss point 
are given in Ref. [16]. Equation (42) may be simplified 


somewhat by combining the summations and weights 


2 
In =D Wy 9(Ey ony) (47) 
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where k = dhey 


We - W x W 


For line integrations, the Gaussian quadrature formula in- 


volves only one summation 


pee NEUTRONIC FIELD EQUATIONS 
The discretization of the spatial domain by the finite 
element method is accomplished by applying equation (37) to 


the governing field equations, using 


o, = <n >® (¥yegjh con es ee (49) 


1. Fuel Region 


The governing equation for the fuel region, equation 
(20) after applying equation (37) becomes 
N. 36 36 
i F ib F 
arf fa = rdrdz-2n f f Noe sp (rDe ao 
eZ ez 


3 Adi ’ oe 
+ a5 (De aoe) +2 poe (1-8) [kZto-bin(T-/Te) J 


= oon wa = 1 — = 
+ iBrap fe A, (t 2 (k2tp)opdt'+iC%e AU) ndrdz = 0 


(50) 
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The second order terms in equation (50) may be re- 
duced to a first order term by application of Green's Theorem 
or equivalently integration-by-parts (See Appendix B). 
Dividing through by 2m and reducing the second order terms 


yields 


ON. 96, IN. O65 , 
yf {else set a ae ANG Pap op (1-8) [kgto-bin(Te/F eo) I] 
le a4 


t 
. a 
-N.IBE ‘| ae ) (kt) o-dt'-N.1C°e** § rdrdz = 0 (51) 
1 aF co F i 
0 


From continuity and boundary conditions the line in- 
tegrals are zero. Now using the approximate functions of 
equation (49) and noting that {pete is not a function of 


space, equation (51) may be written as 
1 ff, ey .e . eft 2 e e 5 2 
v3 {N. Wee rapa tel a DELIN; y) <No alias se Jrdrdz{yp} 


to 
° fe) 2 e el x -r(t-t') ° 1 


e sno At o 
Mfo,r-rprerdotip - AC°e Lf 0%, pravas = 0 (52) 
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The integrations may not be easily carried out with 
a shift to local coordinates and with use of the previously 
derived transformations. Rearranging equation (52) and 


assuming the properties are constant for each time step gives 


De A : C{N. 


* * * e 
Saran 3<dy1 Ns gP*N Jag }<dgg Ny pTrdetLdldedn typ) 


J§ 


ert 
Rap OI-B) (ket) £ L (Ny}<N> rdet[J ]dedn{yp}" 


] Te 
= Ge 8)b- f ‘i In ze = {N, ish >rdet[J]dédn {y-}° 
aF Te F 


lien ‘swe! 
-TBr¢ f(t) i f {N.}<N, pordetLa]dedniy,}* ~icee st 5 i {N, }rdet[J]dEdn 
] ] it . e 
+— sf f {N; }<N.>rdet[d ]dédn {be} = 0 (53) 


es 


Since the element chosen has eight degrees of freedom 
(nodal points), the discretized matrices which result from 
the integration of equation (53) will be 8x8 matrices and the 
forcing function will be an 8x1 vector at the element level. 


Defining the matrices as 


] ] * * * o5 
E a C{N; 2Jqq }<Jqq N; c>tiNs aI50 }<Jo5 N; drdet(d] dEdn 
2 
See «2 
met 2s texe ae PON tS oa 


* 
Gee KON iom ke 22 kk 
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= OS 


8 
where pee WWE? and erates > ri(Ni)y 
i=] 
Tt ye me 
me f {N. I<N, >rdet[J]d&dn = [H3 ss Jexg = a (Ni) K (Ns) rylly (55) 
k=] 
2 
11 H 
Pee iprdetid dein {Ftp = 7 at ) Ply (56) 


ik 
cA | In(Te/Tp) {N, }<N >rdet[d ]dédn = [H4; Joye 
2 
Sees 
aaa aus ice KONG) (NS ) PEW (57) 


To carry out the summation of equation (57)5- the 
temperature Te must be known; however, the temperature is 
exactly what is being sought. To alleviate this problem, a 
linearization is used. 

In the solution technique, the temperature is pre- 
dicted for the next time step. It is this temperature which 
is used for the determination of matrix H4, 


Equation (53) simplifies to 


{DeLH12, ,J-[2,-(1-8) (koto) +28, F(t) ILH3, ;]+25-(1-8)b(H4, 5 T}{v_}° 


Tao Ntre Cwemal soe 
aACTe LE For its elie = 0 (58) 
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The function f(t) is evaluated by summing the 
Values at each time step using the trapezoid rule for numeri- 


cal integration. 


* Ewiheby 
M(t) =oe0 fe ker eiat (59) 
0) 
Defining 
Per 
I,[9(t)] = 5 hjta(ty)+ 9(t;_,)} (60) 
g(t) = e**(k2+p) 
we - time step taken 
Fee 


The function may be expressed as 
5 
Lomeat | 
f(t) =e 2X Teg (ad (61) 
T= 


5 


number of time steps 


es Clad Region 


Following the same procedure as with the fuel equa- 


tion, the discretized form of equation (11) becomes 
ey ol * 1e _ 
{D.[H12; 5] ~ Za lH3; jJ}{v} + “gH3; 510.3 = 0 (62) 


3. Coolant Region 


The governing equation for the coolant region, equa- 


tion (12), may be discretized into the form 


e, | : e_ 
{D. LH12; 5] + Bacolt3, 5 J} beg} + seal eg: = 0 (63) 


ae 


Be 


Once the governing equations have been discretized 
at the element level, they are combined into a system of 
equations at the global level. On the global level the 
governing equation for neutron transport takes the form of 


ae an a CPlixnt} us Eg pe (64) 


nx] 
where [H], [P], and [F] represent the system matrices and 
n is the number of nodal points used in the discretization. 
There are, then, n simultaneous ordinary differential equa- 


tions used to describe the neutron transport problem. 


PeeeneAT TRANSPORT FIELD EQUATIONS 

The spatial domain for the heat transport problem is 
discretized in the same manner as the domain for the neutronics 
problem. The same element matrices previously defined are 


Valid. Let 


ig e e j= 
Te Sige Longe jae eee ae (65) 


1. Fuel Region 


The governing field equation for the fuel, equation 
wee), is discretized by applying equation (37). Using an 
integration by part to lower the order of the second order 


terms allows equation (24) to be written as 


aT a aN, aT, aN, aTe 
FE [Nike 3Z- aia rdrdzt f (Nike ae] dz-S S {kelge- ae + az az] 
7A (0) Veer 
aT 
eZ eNioe + N. (PFC oF Ea an rdrdz = 0 (66) 
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From continuity considerations the line integrals 
are zero except along the boundaries of a region. It is 
assumed that no heat is transferred from the cell in the 
axial direction (boundary condition 5); therefore, the 
first line integral of equation (66) is zero. In the neu- 
tron diffusion problem there was continuity of flux at the 
interfaces so that the line integrals were zero; however, 
the heat transfer at the interfaces is affected by the gap 
and film conductances. The fuel-clad interface condition, 
equation (29), may be rewritten as 


oT - aT 
or... Heap! UF cle? se “Ke ane (67) 


“ke 5 


Substituting equation (67) into (66), dividing by minus one 
and utilizing equation (65) yields 


e e,a =) e e,a e* 
pire eae oes ee - eee ie <M? | dzue! 


g 0 


res kel iN: fey. +{N, 
rz 3 


e e e 
‘ oe S: Stle Jrdrdz{t,} 


> 


e e e 
=u J ele iN. } ie rdrdz{p-} 


fe f f OFC yg (N, }°<N.>Srdrdz (ae = 0 (68) 


e* - element across the interface 
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Transforming to local coordinates and integrating 
by Gaussian quadrature yields the same element matrices as 
given for the neutron flux, equations (54) and (55), except 
for the line integrals between regions. It should be noted 
that the line integrals exist only on the interfaces, along 
which there is a discontinuity of temperature. For the 
fuel equation the interface corresponds to the local co- 
Ordinate €=1. Define 

1 


f {N. I<. >dn = 


Le 
i Ks jdexg =F Dy (Na) (NG) RM (69) 


2 
m 
k=] 
where the N's are evaluated at €=1. Many of the terms of 
Kl will be zero since only the nodes on the &=1 boundary 
Will have shape functions which are non-zero. It is through 


Kl that the temperatures for each region are coupled to- 


gether. Equation (68) may now be written as 


ie Hap (LK Jit, } -(K1]it, he “+k CHI2Z] {723° 


! 
l=) 


~elep(H3]{vp}° + oC -[HB]{t_} = (70) 

In obtaining this equation, it was assumed that material pro- 
perties for each nodal point were constant at each time step. 
Perhaps a better assumption would have been to assume an 
average value for the properties of each element. The dif- 
ference should not be significant, and the assumed constant 


nodal properties were numerically more tractable. 
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2, Clad Region 
Applying equation (37) to the governing field equa- 
tion for the clad, equation (25) gives the discretized form 
of the equation. Assuming no heat transfer in the axial 
direction on the boundaries (Boundary condition 5), equation 


(25) becomes 


ar oN. ae oN. oT. 
- 4 [N. irk. aa dz + c, k elssr + Ee rdrdz 
at 
ype Po Coc N. a rdrdz = 0 (71) 


In the cladding, there are two interfaces along which 
the line integral of equation (71) is not zero, along the 
fuel-clad interface and along the clad-coolant interface. 

For the fuel-clad interface, equation (61) applies. For the 
clad-coolant interface, the interface condition, equation 


(30), may be rewritten as 


ca e Ae oe 
k. —— = hsurf Ole A aay) k (72) 


Along the fuel-clad interface, the local coordinate 


corresponds to &=-1. Define the new element matrix 


em® 


sl ae 
2 


4 ; {N. ic <Ns > “dn = [K2; 5 Iexg 
. k=] 


(Ne), W 


where the N's are evaluated at €=-1. As with Kl, K2 will 
have many zero values because the shape functions are 


evaluated at &=-1. 
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Along the clad-coolant interface, the local coordi- 
Mare,corresponds to €=1] and the Kl matrix is appropriate. 


Substituting equations (67) and (72), equation (71) 


becomes 
thea ea ToT co) I 82 ¥ Pale (Ura TB dz 
oN. aT. oN. aT. 
+ a fk si + Sep rdrdz 
aq. 
ee Pc Coc N. a rdrdz = 0 (Crs) 


The governing equation for the clad region may now be written 
as 


{Dk2}{c }°-[k21itp}°} + rh, peflKI]t © -[KI{t.o}° } 


Pi gap Nur 


+ kTHI2Z]{tp)° + 0, Cy [H3]{t,} = 0 | (74) 


The line integrals of equation (73) affect only nodes which 
are on one of the boundaries; therefore, the nodal inputs 
into Kl and K2 are zero unless the node is on one of the 
boundaries. 
3. Coolant Region 
The field equation governing the coolant may be 
discretized in the same manner as above. After applying 


the Galerkin method and performing an integration by parts 
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on the second order terms, equation (25) 


becomes 
ON. OT N. 9d 
1 co 1 co 
-f{rN. BK eG a eel ct ees ae k lap Se theaige | rdrdz 
Z Ye" Z4 
oT 
co 
Por ifig if Meg? é6 Coco LE rdrdz 
Lie ZA 
al 56 
a8 ae Poo Coco N. la rdrdz = 0 (75) 


The line integral, when evaluated at c, is zero (boundary 
condition 4). When evaluated at b, or correspondingly at 
E=-1, equation (72) is valid and K2 matrix is appropriate. 
All the terms of equation (75) have been defined except the 


flow term. Define 


_ 
e 
4 ef {N, Le <N, na r det[J]d&édn = [H5;s Joy 
Af me 
ate GHaoeTeN, Nene Rewp (76) 
= 


Transforming to local coordinates and integrating reduces 


equation (75) to 


Mb Meee tLkeliz.., } “-[K2]{t, = a + kK. ue iit 


+ tene et Coco ltS]it o} see Creo lt3}{t K9} = 10 C77) 
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Now that the governing equations have been defined 
for each region on the element level, equations (70), (74), 
and (77), they may be assembled into a system equation on 
the global level. The equation will be in the general form 


of 


Wee, ft), * IMI ty 


nxn n x1 2 [GJ iy,tth =u 


n nx] 


Pee CRETIZATION. OF THE, SPATIAL. DOMAIN 

Prior to the numerical solution of the governing equa- 
tions, equations (64) and (78), the spatial domain must be 
divided into a number of elements. For this work the domain 
was subdivided as shown in figure 5. 

Since there is a discontinuity of temperatures at the 
interfaces, as described by equations (29) and (30), a 
noval application of the FEM method was necessary. The com- 
mon practice for handling these "flux' type boundary condi- 
muons 1S to define.a constant reference temperature, T, ; 
as when working with a convection heat transfer problem [17], 
Or to define a known function, as when working with a 
fracture mechanics problem [18]. In either case the refer- 
ence condition was known. The novel application here lies 
in the use of a different field equation to describe the 
reference temperature, e.g., the clad equation (74) is the 
reference condition for heat transfer from the fuel across 


the gap interface. 
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The discontinuity of temperature at the interface neces- 
Sitated another novel application of the FEM. Since there 
is a temperature drop along each interface, a single node 
there is not adequate. In the discretization of the domain, 
two nodes were used for each interface point (for example, 
points 62 and 63 in figure 6). This allows the temperature 
drop due to the gap and film conductances to be taken into 
consideration. Since two nodes are used, the governing equa- 
tions for each region are not directly coupled together. 

The coupling of the regions is accomplished by the "flux" 
boundary or interface conditions since it is assumed that 
any heat flux leaving a region enters the adjacent region. 


Consider a typical set of elements on an interface 


1] 12 13. 14 15 16 


6L 


637l94l45 
et Mies x ¢ 


CLAD—-COOLANT INTERFACE 
FUEL-CLAD INTERFACE 


Figure 6. Finite Element Discretization 
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The coupling terms K1] and K2 may be combined into a system 
K matrix which shows the coupling. The K matrix for the 


Simple set shown is 
iS AS 6 758) 8 ee 18 1A BIS 6 


— 


zz a -a ¢ —¢ c -¢ 

SINC Le eG tw ¢ ease . a= 16/15 
4 b= 4/15 
5 ce = 2/15 
6 d= 1/15 
74 ec -c¢ b -b d -d 

sj —-¢ ¢ -b b =diasd 

13 Ci =¢ d -d b -b 

ma) -¢€ 6¢ -dd -b b 

15 

16 


As can be seen, the nodes on the interface are coupled to 
the adjacent element interface nodes. For example, node 2 


in element I is coupled to nodes 3, 8, and 14 in element II. 


E. OPTIMUM COMPACTING SCHEME 

The system matrices (H, P , etc.) are nxn matrices, 
where n is the number of nodal points used in the discre- 
tization of the domain. In terms of computer storage, these 
matrices may become excessively large if they are stored as 
nxn. There are several techniques available to reduce this 
storage requirement. The most common method is the banded 
storage scheme, whereby only the banded portion of the 


matrices are stored. With judicious numbering of the nodes, 
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considerable savings may be realized. However, it is not 
the optimum storage scheme [8]. 

Since the shape functions, Nis for the poh nodal equa- 
tion are nonzero over only the element containing k, the 
System matrices are not only banded but sparse as well. 

The sparseness is due to the non-consecutive numbering of 
the nodes surrounding the fecal node. The optimum compact 
storage (OCS) scheme compacts the matrices by storing only 
the non-zero elements of the matrices. The implementation 
of the OCS scheme requires two additional integer arrays, 
say JA and NAME. The NAME array identifies the nodal points 
which contribute to each nodal equation. The JA array acts 
as a pointer to indicate where the nodal equation starts in 


NAME. Consider the following simple 2x2 system with nodes 


as indicated 


The NAME array starts with node 1 and identifies the nodes 
mamen contribute to node 1 (i.e., 5, 4, and 2). The NAHE 
array would, then, give the nodes contributing to node 2 
and So forth, so that 


MPQES ene 516)7, O90) eek 1S 20s aS es selene ce JC 


Sess =) <1o42 : 254169 23 65 2 « c.« «-9856> 
and Legit 9 
Shee ea ech 5 4] 4. > 
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The algorithm to assemble the element matrices into a compact 
storage vector is straightforward and represents a significant 
Savings in computer storage [8]. The system matrices are 
stored as a vector rather than a two-dimensional array. For 
example, the value which would be stored in position (1,5) 


Guetne nxn.array is stored in position 2 of the system vector. 
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VI. NUMERICAL SOLUTION 


This section contains a brief description of possible 
solution techniques in addition to the solution technique 
Chosen. Computer subroutines necessary to implement the 


technique are also described. 


Meee SELECTION OF METHOD 

The numerical solution of the system of implicit ordi- 
nary differential equations, equations (64) and (78), may 
be accomplished by any of a number of different techniques 
such as Houbolt's method, Crank-Nicolson's method, Gear's 
method, or implicit Gear's method. It was not the objective 
of this analysis to determine which of the numerical solu- 
tion schemes is the most efficient. Each method has its 
advantages and disadvantages. The Crank-Nicolson method 
is a single-step, implicit equation solver and, therefore, 
does not require storage of previous time solutions. When 
analyzing neutronic problems, the system of equations which 
praises, 1S.commonly very stiff (i.e., a rapid change in flux 
Over a short period of time). The Crank-Nicolson method 
has, in a past work [8], demonstrated difficulty in tracking 
these stiff systems. Gear's method was specifically 
developed for stiff systems and can handle the problem very 
well. However, Gear's method is a multi-step, predictor- 
corrector method requiring storage of previous time solutions. 


In addition to this disadvantage, Gear's method requires the 
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transformation of the developed implicit 0.D.E.'s into an 
eepnert system of 0.D.E.'s.. After this transformation is 
done, the system matrices are no longer sparse or banded, 
thus eliminating the use of the optimum compacting scheme. 
In an effort to overcome these difficulties, Gear's method 
was modified, Ref. [9], to treat the implicit system of 
equations as well as to allow use of the optimum compacting 
scheme. A previous work, Ref. [8], has shown that the 
implicit Gear's method is particularly attractive in solv- 
ing the type problem developed in this analysis. Therefore, 
the implicit Gear's method is used for the solution of the 
system of O0.D.E."s arising -in this analysis. 

No attempt will be made here to give the mathematics in- 
volved in developing the implicit Gear's method. Reference 
[9] may be consulted if details are desired. A listing of 
the computer program developed will be given in the Computer 
Program section. In order to utilize the implicit Gear's 
method, several user supplied subroutines must be developed: 
i) YOLFFUN, 2) JACMAT, and 3), NUITSL. 

B. USER SUPPLIED SUBROUTINES TO IMPLEMENT 
mHE IMPLICIT GEAR'S METHOD 
ie DIFFUN 
Subroutine DIFFUN evaluates equations (64) and (78) 
for a given time and for given values of jy, b, t and t 
Since at each nodal point, i, there is a solution for the 
flux and for the temperature, the solution was set equal to 


DYI and DYII, respectively. In addition to having flux and 
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ey 


temperature at each nodal point, there are also three differ- 
ent regions in the domain which have different governing 
equations. An-integer array, ITYPE, was developed to indi- 
cate for each nodal point whether it was: 0) a fuel node 

not in an interface element, 1) a fuel node in an interface 
element, 2) a cladding node or, 3) a coolant node. Using 
ITYPE, the computer program is directed to a different sec- 
tion depending upon the type of node being considered. After 
all the nodes have been considered, boundary conditions are 
established by changing DYI and DYII for the appropriate 
boundary nodes. Since in this analysis there is continuity 
of flux at the interfaces, special considerations must be 
given to these nodes. At the fuel-clad interface, the value 
of DYI for the clad node was set to the value of the flux 

at that node minus the value of the flux at the adjacent 


mode {(7.e., DYI. = w 


; - V4). Similarly, at the clad-cool- 


F 
ant. interface, the value of DYI for the coolant node was 

set to the value of the flux at that node minus the value 
of the flux at the adjacent node. During the solution of 
the problem, DYI is driven toward zero, which in the limit 


Aner. This is the desired continuity 


forces v; to equal w 
result. 
2. JACMAT 
Subroutine JACMAT, evaluates the Jacobian matrix 
(for Gear's method) at the given time and for the current 
values of the dependent variables. The Jacobian for an equa- 


tion of the type, 
2a ee een (79) 
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may be represented as [19], 


a Sei ey 


where Oo and Ei are coefficients from Gear's method and 
h is the time step. Using the notation of DIFFUN, let DYI 
and DYII represent equations (64) and (78), respectively. 
The Jacobian matrix may, then, be written as (J is called 


PW in JACMAT.) 


aDyI- Yo spi aDNIt Ao aDYII, 
oy Bh ay wees 


(81) 
0 OT 
It is the form of equation (81) which is programmed in 
JACMAT. As in DIFFUN, the integer array ITYPE is used to 
indicate the appropriate section of the program to be 
utilized. The problem boundary conditions must also be 
accounted for in JACMAT. In DIFFUN, the value of DYI or 
DYII was set to zero for constant boundary conditions (i.e., 
zero). This cannot be done in JACMAT since a division by 
zero would occur. For a constant boundary condition at the 
ith node, the value of PW is set to one for the diagonal 
term and zero for all other terms of the ith equation. 
se NUILTSL 

Subroutine NUITSL solves the system of equations 
for the quasi-Newton iterates. In this analysis the system 
is solved using a successive over-relaxation (SOR) method. 
In this work, the optimum amount of over-relaxation was not 


determined. Since no effort was made to find the optimum, 
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ijt was felt a small over-relaxation would be best. The over- 
relaxation factor of 0.02 was used. For small values of 
this factor, the SOR method approaches the Gauss-Siedel 


jteration technique. 
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Vil. PROCEDURE 


In this section, the method utilized to obtain a solu- 
tion is described. The input data necessary to run the 
developed computer program will be documented. 

Prior to initiating a transient overpower excursion, the 
Steady-state conditions for the fuel cell must be known. 
Since the systemof equations which were developed are not 
specifically designed to obtain a steady-state solution, the 
initial steady-state conditions must be part of the input 
data. The initial temperature distribution was obtained 
from the steady-state conditions given in Ref. [7]. 

The axial temperature distribution for the fuel center- 
line, fuel surface, clad, and coolant have been determined 
for several different fuel life cycles [7]. For this analy- 
Sis, the beginning of life cycle for channel 10 was used. 
Although this distribution is somewhat artificial, it should 
be adequate for this analysis. It is the trends of the 
results which are considered important. The distribution 
within the fuel radially is taken to vary as the square of 


the radial distance; then 


re re 
T-(r,z,0) = T-(0,z,0) (1- 2) + Melanesia) 


Within the cladding and the coolant, the initial radial tem- 


Pe~ature distribution is assumed to be constant. 
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The initial flux distribution is assumed to be radially 
Sanstant, a flat flux assumption. In the axial direction 
the flux is assumed to vary as the shape of the sine func- 
tion. The maximum flux, the flux at the axial center, is 
an input parameter. For this analysis, the maximum initial 


iS neutron/em-sec. 


flux was taken to be 10 
To obtain a steady state flux distribution, the value 
of fission cross section for the fuel is varied. A trial- 
and-error method is used until a critical fission cross sec- 
tion, Boe » which gives a steady flux is obtained. 
Once the steady-state conditions have been determined, 


the excess reactivity may be inserted. This starts the 


transient overpower excursion. 


A. INPUT DATA 

The mirst data card contains: the order of Gauss quadra- 
ture, the number of radial elements in the fuel, the number 
of axial elements, number of nodal points in the radial 
G@imection, and the height of the fuel rod. The next cards, 
One for each radial nodal point, contain the nodal radial 
distances. The next cards contain the fuel centerline, fuel 
Surface, clad and coolant temperatures. There is one card 
for each axial node. The next card contains the maximum 
flux. The next four cards contain the physical parameters 
Visted in Table I. The last input cards contain the time, 
end time, estimated initial time step, minimum time step, 


and maximum time step. A sample data deck is shown in figure 


is 
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TABLE I. 


Fuel Diffusion Coefficient (DCF) 
Doppler constant (B) 

Energy release per fission (E) 
Fraction of delayed neutrons (BETA) 


Fraction of delayed neutrons for 
the ith group (BETAI) 


Decay constant for ith delayed neutron 
group (DCLAMI) 


Initial flux for delayed neutrons 


Average number of neutrons released 
per fission (ANU) 


Neutron velocity (VEL) 
Fuel absorption cross section (SIGAF) 


Critical fuel fission cross section 
(SIGFF) 


Blanket absorption cross section (SIGAB) 
Blanket fission cross section (SIGFB) 
Step reactivity input (RHOA) 

Ramp reactivity input (RHOB) 

Fuel density (DENF) 

Clad diffusion coefficient (DCC) 

Clad specific heat (CPC) 

Clad density (DENC) 

Clad thermal conductivity (TKC) 

Clad absorption cross section (SIGAC) 
Coolant diffusion coefficient (DCCO) 
Coolant absorption cross section (SIGACO) 
Coolant flow velocity (VCO) 

Surface heat transfer coefficient(HSURF) 
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Physical Parameters 


0.93 [cm) 
0.006 
7.652x10 
0.0064 


-12 


0.0064 


0.0784 
1x19 !9 [neutron/cm¢sec] 


2.44 

4.8 x10° [cm/sec ] 
0.088 [em7!] 
0.0586875 [cm™'] 


0.0 Com7!] 

0.0 [em !] 
Variable 

Variable 

10.9 [gm/cm?] 

1 Lem] 

0.12 [cal/gm °C] 
8.0 [gm/cm?] 

0.0526 [cal/cm sec °C] 
0.0015 [em™'] 

1.55 [cm] 

0.00004 [cm™!] 
396.0 [cm/sec ] 

0.7 [cal/cm*sec oC] 


[cal/fissions ] 
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VOT RESUS 


When using the finite element method, one of the first 
considerations must be given to the convergence of the 
method. To determine convergence, the results for a given 
point are compared for different finite element discretiza- 
tions. As shown in figure 8, the results are comparable 
but no definite claim of convergence can be made. However, 
for this work, it was felt that these results were adequate. 
It was not the object of this analysis to arrive at the 
"final" result; it was the trends and methods that were of 
interest. Since the 66-element mesh appears to give a fair 
approximation of the results, the 66-element mesh was used 
as the discretized domain. 

The next item of consideration was the determination of 
a neutronic steady-state condition. This proved to be a 
werys time consuming task. The fission cross section for the 
fuel was varied by a trial-and-error method in an attempt 
to find the critical cross section which would give a steady 
State. As may be seen in figure 9, a change in cross sec- 
tion of less than one percent significantly affected the 
State of the problem. It was felt that the critical value 
Was between the values of 0.05875 and 0.058625. Time did 
not permit investigation for critical value; therefore, 
it was assumed that the value for the critical fission cross 

cr 


Section was half way between the values (i.e., Le = 0.0586875) . 
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Even if this value is in error, which it most probably is, 
the net effect would only be a small decrease or increase 
in the proposed reactivity insertion. The reactivity inser- 
tion would, then, be an approximation of the actual reactiv- 
ity of the problem. 

The first test problem considered was a step increase 
in reactivity of approximately ten dollars. For the uranium 
dioxide fuel, one dollar of reactivity was taken to be 
meooes. Figure 10 shows the time history of the flux at 
the center of the fuel rod. Figure 11 gives the correspond- 
ing temperature profile. The temperatures were taken at the 
hottest point of each region at the axial center (i.e., the 
fuel centerline, clad inside surface, and coolant inside 
surface). As seen on the fuel temperature time history, the 
fuel rapidly reaches the fuel melting point. The model 
developed does not take into consideration melting of the 
fuel. This melting would tend to decrease the effect of the 
transient. The problem was allowed to continue despite this 
inconsistency in the mathematical model. A short time after 
fuel melting, the inside surface of the clad reaches its 
melting point. The temperature in the coolant experienced 
What is felt to be a numerical phenomenon. The coolant tem- 
perature decreased prior to the small rise at the end of the 
transient. Intuitively, this decrease does not seem to be 
realistic. A similar occurrence was observed while conduct- 
ing sample tests on the developed computer program. Refer- 


ence [20] reported the same phenomenon. It is felt that this 
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phenomenon is a quirk of the finite element method. The 
radial and axial distributions of the neutron flux are pre- 
sented in figures 12 and 13 for time equal to 7.39 x Tee 
Seconds. The distributions are, basically, as anticipated. 
The neutron flux peaks slightly before the axial center. 

It was expected to peak at the axial center. The radial 

and axial temperature distributions for the same time are 
presented in figures 14 and 15. As with the flux, the tem- 
perature profiles were, basically, as expected. The fuel 
temperature peaks slightly below the expected location, 

most likely in response to the peak in axial fluxes. The 
coolant unexpectedly drops near the outlet of the fuel rod. 
The finite element method characteristically has some prob- 
lems on the boundaries of the domain; this may account for 
the drop in coolant temperature. 

The temperature of the fuel and cladding do not appear 
to be as closely coupled as anticipated. As seen in figure 
14, a significant increase in fuel temperature has resulted 
in a relatively small clad temperature increase. As noted 
in figure 11, there appears to be a time lag in temperature 
response for each region which may account for part of the 
apparent temperature disagreement. It is felt that the 
temperatures should be more closely related, which indicates 
a higher gap heat transfer coefficient should be utilized. 
Since values of the gap heat transfer coefficient were 
assumed, it is not unreasonable to believe the values used 


are too low. 
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Time did not permit investigation of other reactivity 
insertions. Other reactivity inputs may be investigated by 
students in the future. 

This work does not represent a solution to the very com- 
plicated nuclear reactor problem. It does represent an 
application of a numerical technique which is relatively new 
to nuclear applications. Methods for implementing the finite 
element method have been discussed, and a computer code has 


been developed for the simplistic model considered. 
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IX. RECOMMENDATIONS 


For the model developed, perhaps the most important item 
Eoepursue is the critical fission cross section. A better 
determination of this value is necessary so that the reactiv- 
ity insertion is more accurately known. Different test cases 
mor the prompt critical and prompt subcritical reactor could 
then be conducted. 

In further developing the model, more consideration 
Should be given the gap heat transfer coefficient. As noted 
in the results, the value used appears to be too small. 
Sample problems for different gap heat transfer coefficients 
would give a better indication of the values to use. 

Melting of the fuel during the transient would probably 
be the next major improvement on the model. With relatively 
few changes, the model could be adapted to allow melting 
element by element. This, too, would be an approximation 
but, still, an improvement to the model. Perhaps at the 
Same time, a simplified model to take into consideration the 
fuel restructuring could be implemented. 

Another improvement would be to consider reactivity 
feedbacks in addition to the Doppler feedback. Sodium void- 
ing and fuel rod expansion are two of the more important 
feedback effects to consider. 

On the numerical side, probably the most important thing 


to do would be to run the computer program on the 
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"H-compiler", which optimizes the program. However, on 
several runs using the H-compiler, erroneous results were 
Obtained. With sufficient time, this could be corrected 

to allow use of the H-compiler. The use of the H-compiler 
results in a savings in computer time. The present program 
runs on a "G-compiler" and takes excessive amounts of com- 
puter time (two to four hours per run). 

In addition to this, the optimum over-relaxation factor 
in the implicit Gear's method could be determined by trial- 
and-error. 

Implementation of these recommendations should enhance 


the analysis and lead to a more efficient computer code. 
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APPENDIX A 
DEVELOPMENT OF TRANSFORMATIONS 


The Jacobian matrix [J] may be written (8) for two 


dimensions as 


N N 
2 ee X Tee 
1= = 
[o> ‘ j (Al) 
du eee yy a 1 


= 1 d- =b 
Se det[A] Bc 4 


Applying this fact to equation (Al) gives 


N N 
2a NaAh2, > a ees 
i=] 


=] 
-1 
idi~ = TCR N N (A2) 
nee ueml treet daha 
* * 
— | 44 di 
* * 
Joy Joo 
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From matrix algebra 


N 
det [J] >> Satz Mo rR Ny eee AS) 


i=] i=] i=] 


The derivatives of the shape functions may be found from 


equations (38). 


Ny og = g(ltn) (2E4n) Nyon = Ge lte) (2nt6) 

Nog = -6(14n) No, = gle") 

Nag = q(ltn) (2é-n) Nao = gle) (2n-8) 

Nae to glen’) Ng, 7-1-8) (A4) 
Neg = q(I-n) (26+n) Ne = ll-8) (2nt6) 

Ne eg = 7E(1-n) Ne = > p18") 

Nog = p(Ion) (26-n) Noon = q1t8) (2n-6) 

Neg = ellen’) Ng, = -n(14é) 


If one now considers an arbitrary element 
z 


PAY) 


6 
| 
| | 
| l 
| 
| 
| | 
r r 


| 
| 
| 
| 
| 
| 2 3 r 
with the midside nodes exactly at the midpoint (not a neces- 


Sary criteria for the FEM) and substitutes into equation (A3), 
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the result is 


aah (iss 
a ea Te ve 


Substituting into equation (A2) and using (A5) will yield 


J ee liz =Z2- )) ="2/%s=0 (A6) 
| fee So ae ee 
* 
ie ee (A7) 
aaa 
Joy) Os (A8) 
and fed as 2 le er.) ) = 2725-2 (A9) 
22 fee oe cen lah 


The inverse of the Jacobian matrix now becomes 


2/r3-r, 0 
[J] = (A10) 


0 2/Z3-Z, 


For integration along a line, the transformation used 


1S 
dz = det[J']dn (A111) 
In this case 
N 
det[J'] = 2 Hehe (A12) 
1= 


Again considering the arbitrary element and substituting 


(A4) into (A12) will result in 


TAD 7h e 
det[J'] = = ie > (A13) 


APPENDIX B 
REDUCTION OF SECOND ORDER TERM 


The second order term of the governing field equations 
may pe reduced to first corder by integration by parts. 


Consider, for example, 
-- = 
Sf «a NL = = ( rD =a sy NO =e rdrdz (B1) 


which may be expanded as 


2, 
Pf 070 2 + nye 0 Mo nro AG + nye 20 3 drdz 
9r° 
(B2) 


Integrating just the second order terms by parts will yield 


2 
acy i ay 
Sf cee oraraz =) INjrD srl p42 
If 74 of Zé 


ee Semon, + YN, = + rD jal drdz  (B3) 


and 


Z F 

Cm e aD 
Sf v,0 as rdrdz = fnj0 3 ey rdr Sf ao, Oty asl rdrdz_ (B4) 
2 r 


Substituting the results of (B3) and (B4) into equation (B2) 


will give 


F aN. 
fo rD cles + [N, D oy rd ff = at = A rdrdz_ (B5) 
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APPENDIX C 
LIST OF RELATIONS FOR MATERIAL THERMAL PROPERTIES 


eee FUEL (U0, ) 


me Specific Heat, Ref. [19] 
Ope [18.45#2.431x107°T=2. 272x107 °1°]/270.07 
feal/gm °C] 


i; =<, aE 
Ze Thermal Conductivity, Ref. [5] 


Tke = [1-2.5( 1-75) Ixlyqeee + 4.79x10 


14303 
3541 T 


1x0.235 
feal/em “sec 7 C] 
Prp - percent theoretical density 


b - °K 


B. CLAD (Stainless Steel) 


Properties are assumed to be temperature independent, 
and average values from Ref. [5] were used for the 
clad properties. 


C. COOLANT (Liquid Sodium) 


ie ‘Specific Heat, Ref. [5] 


4 Dae 


Cy = 0.34574-0.79226x10 “T+0.34086x10 ‘T 
co 


feal/gm-° C4 
fear 


ee Density, Ref. [5] 


3 2 


T-0. 2872x107 °T 


Peo. [59.566-7.9504x10° 


-953 


Cc 


+ 0.06035x10 J] x 0.01601 [gm/cm?] 


io oF 
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3. Thermal Conductivity, Ref. [5] 


-2742.0914x107°T2] 


3 


to = [54.306-1.878x10 
x 4. 134x110 - Leal/em see °C] 


T= °F 


Mee SURFACE HEAT=TRANSFER COEFFICIENT 


| Lisa De Pearce pe 0.8 
Meare aM [7.0+0.025.( Tk ) ] 
(ol 
[BTU/hr ft* °F] 
Koo =_ LB Ushr Ft. 2F] 
3 
Pos - [ibmé ft] 
as - [ft/sec] 
C - [BTU/1bm °F] 
Poo 


De - equivalent diameter [ft] 
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THE INITIAL TEMPERATURE DISTRIEUTION 
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